/*  scommult.cpp  by K.Tsuru  */
// function ID = 903
/********************************************
SComplex class
multiplication (*this)*z = (a + bi)*(c + di)
********************************************/
#ifndef SN_H
#include "sn.h"
#endif
SComplex& SComplex::operator*=(const SComplex& z) {
    if ( IsZero(903) ) return *this;    // a = b= 0
    if ( z.IsZero(903) ) {              // c = d = 0
        SetZero();  return *this;
    }
    // d = 0, (a+bi)*c
    if (z.im.Sign(903) == SNumber::ZERO) return (*this) *= z.re; // d=0 [operator*=(const SDouble& x) is used]
    // c = 0, (a+bi)*di = -b*d +(a*d)i
    if (z.re.Sign(903) == SNumber::ZERO) { // c=0, (a+bi)*di = -bd+i ad = -im*z.im + i re*z.im
        const SDouble r = -im * z.im;
        im = re * z.im;
        re = r;
        return *this;
    }

//  if(&x == &y) {  // z * z can be speed-up by z = a + ib, z^2 = a^2-b^2 + 2iab
    if(*this == z) {    // Overhead is very small comparing above.
        SDouble temp(0.0);  // (a+ai)^2 = 2a^2i

        if(re != im) temp = re * re - im * im;
        im = re * im;   im = DsMult(im, 2);
        re = temp;
        return *this;
    } else {
    /******************************* 
     (a + bi)*(c + di) = a*c - b*d +(a*d + b*c)i
     = a*(c -d) + d*(a -b)+{b*(c+d) + d*(a -b)}i
     = (re + i*im)*(z.re + i*z.im)
     Let re = a, im = b, z.re = c, z.im = d, temp = d*(a -b).
    ***********************************/
    const SDouble temp = z.im * (re - im); // three times *operator()
        re = re * (z.re - z.im) + temp;  // a*(c-d)+temp
        im = im * (z.re + z.im) + temp;  // b*(c+d)+temp
        return *this;
    }
}